library(tidyverse) #ggplot2, dplyr, tibble, etc.
library(plotly) #use ggplotly() for interactive plots with scoll-over IDs
library(knitr) #use kable() to make formatted tables
library(Rtsne)
library(glmmTMB)
library(lme4)
library(sjPlot)
library(ggpubr) #ggarrange() figure wraps
library(plyr)
library(ggeffects) #ggpredict
library(rstanarm) #stan models
library(shinystan) #stan model evaluation >launch_shinystan_demo()
library(loo) #loo() to compare fits between bayesian models
library(MuMIn) #dredge, model averaging
library(AICcmodavg) #aictab() for AICc model comparisons
library(RColorBrewer)

full mouth color dataset (mc)

mc=tibble(read.csv("CB-mouthColor.csv", h=T))


full nestling morphometrics dataset (nm)

nm=tibble(read.csv("nestlingMorph.csv", h=T))
nm[nm==0] <- NA


Morphometrics dataframe simplified for PCA

#subset id, billD, billW, TEC, skull, tarsus, and age
morphSimple.nm <- nm[,c(6,7,11,16:18,20,23)]
morphSimple.nm <- morphSimple.nm %>% 
  filter(between(age, 23, 27))
morphSimple.nm <- drop_na(morphSimple.nm)
#morphSimple.nm <- morphSimple.nm[-c(640,840,1186,1922,2196),]
morphForPCA <- as.matrix(morphSimple.nm[,c(3:7)])


Size PCA (same metrics as Townsend et al. 2010)

sizePCA <- prcomp(morphForPCA)

sizePCA.plot <- tibble(PC1=sizePCA$x[,1],
                        PC2=sizePCA$x[,2],
                        id=morphSimple.nm$id)
summary(sizePCA)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5
## Standard deviation     4.0329 2.2578 1.7519 1.44284 0.58696
## Proportion of Variance 0.6056 0.1898 0.1143 0.07751 0.01283
## Cumulative Proportion  0.6056 0.7954 0.9097 0.98717 1.00000
sizePCA
## Standard deviations (1, .., p=5):
## [1] 4.0328741 2.2577553 1.7518941 1.4428431 0.5869643
## 
## Rotation (n x k) = (5 x 5):
##               PC1         PC2         PC3        PC4         PC5
## tarsus -0.8008200 -0.08478806  0.15408891 -0.5716853  0.03051037
## billD  -0.1357370  0.04076158 -0.06506740  0.1141950 -0.98114203
## billW  -0.1200268  0.05105806 -0.98367877 -0.1007184  0.07223945
## TEC    -0.4511211  0.68408742  0.04531459  0.5507967  0.15193321
## skull  -0.3497801 -0.72150326 -0.04843046  0.5887395  0.09015093
weightByTarsusPlot <-morphSimple.nm %>% 
  ggplot(aes(x=tarsus, y=weight, label=id))+
  geom_point(shape=1)+
  geom_smooth()
ggplotly(weightByTarsusPlot)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
massTarsus.mdl <- lm(weight~tarsus, data = morphSimple.nm)
plot(massTarsus.mdl)

morphSimple.nm$massTarsusResiduals <- residuals(massTarsus.mdl)
joinResiduals.df <- data.frame(id = morphSimple.nm$id, 
                               residualsMassTarsus = morphSimple.nm$massTarsusResiduals)
new.mc <- left_join(mc, joinResiduals.df, by = "id")
mc$residualsMassTarsus
## Warning: Unknown or uninitialised column: `residualsMassTarsus`.
## NULL

Residuals (mass by tarsus) model for Saturation 1

#checked age * tarsus interaction...not significant
# sat1_resid.mdl <- lmer(sat1 ~ age + residualsMassTarus.x + sex + bodTemp1 + ambTemp1 + timeOut1 + (1|family), 
#                  data = new.mc, REML = FALSE)
# 
# summary(sat1_full.mdl)

Full model for Saturation 1

#checked age * tarsus interaction...not significant
sat1_full.mdl <- lmer(sat1 ~ weight + age + tarsus + sex + bodTemp1 + ambTemp1 + timeOut1 + (1|family), 
                 data = mc, REML = FALSE)

summary(sat1_full.mdl)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: sat1 ~ weight + age + tarsus + sex + bodTemp1 + ambTemp1 + timeOut1 +  
##     (1 | family)
##    Data: mc
## 
##      AIC      BIC   logLik deviance df.resid 
##   1322.8   1356.7   -650.4   1300.8      151 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9805 -0.4147  0.1147  0.5441  1.7181 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  family   (Intercept) 133.65   11.561  
##  Residual              94.32    9.712  
## Number of obs: 162, groups:  family, 88
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 190.47318   55.82130   3.412
## weight        0.09456    0.03766   2.511
## age          -3.51493    0.62139  -5.657
## tarsus       -1.34739    0.64571  -2.087
## sexF          5.74332    9.94909   0.577
## sexM          3.49846    9.74417   0.359
## bodTemp1      1.26644    1.42201   0.891
## ambTemp1      0.34515    0.40305   0.856
## timeOut1      0.01686    0.09879   0.171
## 
## Correlation of Fixed Effects:
##          (Intr) weight age    tarsus sexF   sexM   bdTmp1 ambTm1
## weight    0.447                                                 
## age      -0.157 -0.197                                          
## tarsus   -0.395 -0.748  0.075                                   
## sexF     -0.133 -0.006  0.006  0.121                            
## sexM     -0.157 -0.029  0.004  0.082  0.973                     
## bodTemp1 -0.780 -0.114 -0.123 -0.162 -0.104 -0.046              
## ambTemp1  0.261 -0.065 -0.016  0.209 -0.080 -0.093 -0.496       
## timeOut1 -0.087 -0.028 -0.166  0.076 -0.013 -0.023  0.062 -0.065

Confidence intervals for Full Saturation 1 model

confint(sat1_full.mdl)
## Computing profile confidence intervals ...
##                    2.5 %       97.5 %
## .sig01        8.89525628  14.58041653
## .sigma        8.29391209  11.52047360
## (Intercept)  80.19693221 301.08102026
## weight        0.01768514   0.17124532
## age          -4.75215882  -2.28733611
## tarsus       -2.62435838  -0.07139622
## sexF        -14.06523964  25.62651294
## sexM        -15.78560888  22.84723083
## bodTemp1     -1.58278201   4.08843893
## ambTemp1     -0.48547219   1.15254538
## timeOut1     -0.18884662   0.22891007
#plot model with sex
plot_model(sat1_full.mdl,
           title = "Saturation 1")


#plot model without sex to view the smaller confidence intervals 
plot_model(sat1_full.mdl, 
           title = "Saturation 1",
           rm.terms = c("sex [F]","sex [M]"))


Reduced models: fixed effects removed based on confidence intervals from the full model

#sex random effect 
sat1_rsex.mdl <- lmer(sat1 ~ weight + age + tarsus + bodTemp1 + ambTemp1 + (1|sex), 
                 data = mc, REML = FALSE)
#remove timeOut
sat1_r1.mdl <- lmer(sat1 ~ weight + age + tarsus + sex + bodTemp1 + ambTemp1 + (1|family), 
                 data = mc, REML = FALSE)
#remove ambTemp
sat1_r2.mdl <- lmer(sat1 ~ weight + age + tarsus + sex + bodTemp1 + (1|family), 
                 data = mc, REML = FALSE)
#remove bodTemp
sat1_r3.mdl <- lmer(sat1 ~ weight + age + tarsus + sex + (1|family), 
                 data = mc, REML = FALSE)
sat1_candset <- c(sat1_full.mdl,sat1_rsex.mdl,sat1_r1.mdl,sat1_r2.mdl,sat1_r3.mdl)
mdls <- c("full", "randSex", "r1", "r2", "r3")
aictab(sat1_candset, modnames = mdls, sort = TRUE)
## 
## Model selection based on AICc:
## 
##          K    AICc Delta_AICc AICcWt Cum.Wt      LL
## full    11 1324.52       0.00      1      1 -650.38
## r1      10 1560.88     236.36      0      1 -769.84
## randSex  8 1608.09     283.57      0      1 -795.65
## r2       9 1810.46     485.94      0      1 -895.81
## r3       8 1831.02     506.49      0      1 -907.17

AICc model selection shows that the full model is better fit than any of the reduced models.

In summary, weight is positively associated with saturation 1. Age and tarsus are both negatively associated with saturation 1. Sex, ambient temp, and body temp are not significantly associated with saturation 1.